### texreg method for polr objects
extract.polr <- function(model, include.thresholds = TRUE, include.aic = TRUE,
                         include.bic = TRUE, include.loglik = TRUE,
                         include.deviance = FALSE,
                         include.nobs = TRUE, include.lrtest = TRUE,
                         include.lrtest.divided = TRUE,
                         include.score.test = TRUE, ...) {
  s <- summary(model, ...)
  lrtest.chi <- round(lmtest::lrtest(model)$C[2], 3)
  lrtest.pval <- round(lmtest::lrtest(model)$P[2], 3)

  lrtest.divided.chi <- round(lmtest::lrtest(eval(parse(text = paste(deparse(substitute(model)),".nodivided",sep = ""))))$C[2], 3)
  lrtest.divided.pval <- round(lmtest::lrtest(eval(parse(text = paste(deparse(substitute(model)),".nodivided",sep = ""))))$P[2], 3)
  
  score.mod.po <- eval(parse(text = paste("vglm(",sub("factor", "ordered", as.character(model$call)[2]), ", data = ", as.character(model$call)[3], ", family = 'cumulative'(parallel = T))", sep = "")))
  
  score.mod.nopo <- eval(parse(text = paste("vglm(",sub("factor", "ordered", as.character(model$call)[2]), ", data = ", as.character(model$call)[3], ", family = 'cumulative'(parallel = F))", sep = "")))
  
  score.mod.test <- deviance(score.mod.po) - deviance(score.mod.nopo)
  score.mod.pval <- pchisq(score.mod.test, df = df.residual(score.mod.po) - df.residual(score.mod.nopo), lower.tail = FALSE)
  
  tab <- s$coefficients
  zeta.names <- names(s$zeta)
  beta <- tab[!rownames(tab) %in% zeta.names, ]
  thresh <- tab[rownames(tab) %in% zeta.names, ]
  
  if (sum(!rownames(tab) %in% zeta.names) == 1) {
    beta <- t(beta)
    rownames(beta) <- rownames(tab)[!rownames(tab) %in% zeta.names]
  }
  if (sum(rownames(tab) %in% zeta.names) == 1) {
    thresh <- t(thresh)
    rownames(thresh) <- rownames(tab)[rownames(tab) %in% zeta.names]
  }
  
  threshold.names <- rownames(thresh)
  threshold.coef <- thresh[, 1]
  threshold.se <- thresh[, 2]
  threshold.zval <- thresh[, 1] / thresh[, 2]
  threshold.pval <- 2 * pnorm(abs(threshold.zval), lower.tail = FALSE)
  
  beta.names <- rownames(beta)
  beta.coef <- beta[, 1]
  beta.se <- beta[, 2]
  beta.zval <- beta[, 1] / beta[, 2]
  beta.pval <- 2 * pnorm(abs(beta.zval), lower.tail = FALSE)
  
  if (include.thresholds == TRUE) {
    names <- c(beta.names, threshold.names)
    coef <- c(beta.coef, threshold.coef)
    se <- c(beta.se, threshold.se)
    pval <- c(beta.pval, threshold.pval)
  } else {
    names <- beta.names
    coef <- beta.coef
    se <- beta.se
    pval <- beta.pval
  }
  
  n <- nobs(model)
  lik <- logLik(model)[1]
  aic <- AIC(model)
  bic <- BIC(model)
  dev <- deviance(model)
  lrt <- lrtest.chi
  lrp <- lrtest.pval
  lrt.d <- lrtest.divided.chi
  lrp.d <- lrtest.divided.pval
  
  gof <- numeric()
  gof.names <- character()
  gof.decimal <- logical()
  if (include.aic == TRUE) {
    gof <- c(gof, aic)
    gof.names <- c(gof.names, "AIC")
    gof.decimal <- c(gof.decimal, TRUE)
  }
  if (include.bic == TRUE) {
    gof <- c(gof, bic)
    gof.names <- c(gof.names, "BIC")
    gof.decimal <- c(gof.decimal, TRUE)
  }
  if (include.loglik == TRUE) {
    gof <- c(gof, lik)
    gof.names <- c(gof.names, "Log Likelihood")
    gof.decimal <- c(gof.decimal, TRUE)
  }

    if (include.deviance == TRUE) {
    gof <- c(gof, dev)
    gof.names <- c(gof.names, "Deviance")
    gof.decimal <- c(gof.decimal, TRUE)
  }
  
  if (include.lrtest == TRUE) {
    gof <- c(gof, lrt, lrp)
    gof.names <- c(gof.names, "Likelihood Ratio Test", "%Likelihood Ratio p-value")
    gof.decimal <- c(gof.decimal, TRUE, TRUE)
  }

  if (include.lrtest.divided == TRUE) {
    gof <- c(gof, lrt.d, lrp.d)
    gof.names <- c(gof.names, "Likelihood Ratio Test of Significance of Divided Government", "%Divided Likelihood Ratio p-value")
    gof.decimal <- c(gof.decimal, TRUE, TRUE)
  }

  if (include.score.test == TRUE) {
    gof <- c(gof, score.mod.test, score.mod.pval)
    gof.names <- c(gof.names, "Score Test", "%Score Test p-value")
    gof.decimal <- c(gof.decimal, TRUE, TRUE)
  }
  
  if (include.nobs == TRUE) {
    gof <- c(gof, n)
    gof.names <- c(gof.names, "Number of Observations")
    gof.decimal <- c(gof.decimal, FALSE)
  }
  
  tr <- createTexreg(
    coef.names = names,
    coef = coef,
    se = se,
    pvalues = pval,
    gof.names = gof.names,
    gof = gof,
    gof.decimal = gof.decimal
  )
  return(tr)
}
